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Abstract 

Compressed manifold modes are locally supported analogues of eigenfunc¬ 
tions of the Laplace-Beltrami operator of a manifold. In this paper we describe 
an algorithm for the calculation of modes for discrete manifolds that, in ex¬ 
periments, requires on average 47% fewer iterations and 44% less time than 
the previous algorithm. We show how to naturally order the modes in an 
analogous way to eigenfunctions, that is we define a compressed eigenvalue. 
Furthermore, in contrast to the previous algorithm we permit unlumped mass 
matrices for the operator and we show, unlike the case of eigenfunctions, that 
modes can, in general, be oriented. 


1 Introduction 

The eigenvalues and eigenfunctions of the Laplace-Beltrami operator (LBO) on a 
discrete manifold have found many applications in geometry processing, for exam¬ 
ple, in shape matching, remeshing (such as quadrangulation), smoothing, and shape 
identification, see for example, mi- 

In analogy with a Fourier basis, the eigenfunctions of the LBO form a basis, 
called the manifold harmonic basis (see [II]), of the functions on the manifold. 
One major drawback of this basis is that it does not relate, at least in an intuitive 
way, to observable features of manifolds. In the pioneering work [nnn it is shown 
that one can produce on a set of functions with localized support, i.e., compactly 
supported, which will function as a basis. They named these compressed modes. 
They proposed that these may be used as a natural basis for solving PDFs and 
suggested that they could be extended to general discrete manifolds. This general¬ 
ization was made in m and were called compressed manifold modes (CMMs). 
These CMMs are locally supported and, crucially, they seem to be supported at 
important features of the manifold. For example in the top line of Figure one can 
see that the collection of six CMMs is supported on the arms, legs, head and lower 
torso. That is, features that a human would identify are detected by the modes. 
The six modes pictured on the less symmetric Aquarius the Water Carrier, are also 
supported on regions which a human may identify as regions of interest. Twenty 
modes are given later in Figure for the classic teapot. 

Similarly the support for 15 modes on a L-shaped mesh is local and identifies 
significant features such as corners as shown in Figure 
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Figure 1: Local support and identification of natural features. 





Figure 2: Support on a L-shaped mesh, red is positive and blue is negative. 
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In [To] the authors show that, in contrast to the eigenfunctions and eigenvalues of 
the LBO, the CMMs are robust with respect to noise and holes. Examples are given 
where shape matching is achieved in the presence of noise and partial sampling of 
meshes. Furthermore, the theory in [TT1[T2] is developed as a way of solving partial 
differential equations. Hence, CMMs are potentially important objects of study in 
a diverse range of subjects. 

Two significant problems of compressed manifold modes that are identified in 
m are speed and ordering. They take longer to generate than eigenfunctions and 
no method for ordering was given. Further problems are that the algorithm in m 
is restricted to an LBO with diagonal mass matrix and it is noted that there is 
ambiguity of sign of the modes in the same way that unit eigenfunctions of a matrix 
are only defined up to sign. 

This paper deals with all these problems. The main contributions are as follows: 

(i) . We adapt an algorithm from [5] to speed up the algorithm given in (TOj for 

computation of the modes. This accelerated algorithm is described in Sectionj^ 
and experiments show an average 44% reduction in the time taken with a 47% 
decrease in the number of iterations required. The accuracy of the accelerated 
version is measured in Section [H 

(ii) . We describe a natural ordering for the modes in Section that is entirely 

analogous to the ordering of eigenvalues of a matrix. 

(iii) . The new algorithm allows non-diagonal mass matrices. The explanation of 

this is given in Section 

(iv) . Section]^ shows how, in contrast to the case of eigenvectors, the modes can be 

oriented, i.e., can be given a sign. This can be seen by comparing the colours 
in Figure]^ with that of Figure 13 in m- The modes in Figure]^ have been 
‘flipped’ where necessary. 

The Matlab code of the implementation of the algorithms is freely available 
under a Creative Commons Attribution-NonCommercial-ShareAlike (CC BY-NC¬ 
SA) license. It can be found at http://wwwl.maths.leeds.ac.uk/~khouston/, 
Acknowledgements: My thanks to Thomas Neumann for helpful discussions 
and comments. The meshes of humans come from the SCAPE dataset, [1] and the 
Aquarius mesh is from the EPFL Computer Graphics and Geometry Laboratory. 
The Algorithms Latex bundle by Rogerio Brito was used in preparation of the 
manuscript. 

2 Background, previous and related work 

In recent years there has been much interest in the eigenvalues and eigenfunctions 
of the discrete Laplace-Beltrami operator. For a smooth manifold M the classical 
differential geometry Laplace-Beltrami operator, denoted A, has a set of eigenfunc¬ 
tions {^k} and associated eigenvalues, {A/^}, determined by 

Aifk = Xk^k 

where k e N and Xk G M. See [2]. The self-adjointness of A implies that the 
eigenvalues are real and that the eigenfunctions are orthogonal with respect to the 
L 2 -inner product: (/, g) = fg. 

In the discrete case we denote the LBO by L and decompose it as L = A~^W 
where A is a mass matrix, i.e., a matrix which relates to the area/volume around 
the vertices of the discrete manifold, and IT is a weight matrix, such as the cotan 
matrix, (see HE] for further references). 
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Both A and W are symmetric N xN matrices where N is the number of vertices. 
The matrix L = A~^W is in general not symmetric but is self-adjoint with respect 
to the I/ 2 -inner product given by {v^ w) = v^Aw^ where denotes the transpose 
of the vector v. This implies that its eigenvalues are real and its eigenvectors are 
orthogonal with respect to the inner product. The eigenvalues and eigenvectors of 
L are solutions of the symmetric eigenvalue problem 

Wip = XAip. 

In this paper we write K eigenfunctions as the columns of the N x K matrix 

In [TT] compressed modes are introduced as a solution of a minimization prob¬ 
lem involving the Hamiltonian operator H = — (1/2)A + H(x), where H is a poten¬ 
tial function, with what is called an Li-regularizing term, that is, (l//i)X]j 
for modes 2 pj where /i G M is a parameter controlling how localized the modes are. 

We produce a set arising from the variational problem 

min (— such that = Id 

^ j=i / 

where T is the N x K matrix given by arranging the in columns. 

The generalization, called compressed manifold modes, described in m is 
as follows. For K G N and /i G M, the first K compressed manifold modes are the 
columns of the matrix where ^ is determined by the constrained optimization 
problem 

nnnTr + /i| |^| |i such that A^ = Id . 

Here ||^||i is the sum of the absolute values of entries of (Note that, rather 
confusingly, the compression parameter is l//i in [TT] and /i in [TO] . We will follow 
the latter as our algorithm is a generalization of theirs.) A mode is an analogue of 
the Laplace-Beltrami eigenfunctions. 

The parameter // is a measurement of the compression of the support of the 
modes. A large /r means large compression of support, i.e, the support gets smaller, 
and a small fi has small compression, i.e., large support. 

3 Accelerated ADMM Algorithm 

This section contains the first contribution of the paper - an accelerated version 
of the algorithm presented in m to calculate the compressed manifold modes. In 
the tests this new algorithm required 47% fewer iterations to reach convergence and 
resulted in a 44% time improvement. 

In m, Ozolins et al propose a Splitting Orthogonality Constraint (SOC) al¬ 
gorithm (described in detail in [7]) for calculating compressed modes in the case 
that the manifold is BA. In (TOl, Neumann et al show that the SOC algorithm does 
not function well for more general manifolds and to counteract this they propose 
an Alternating Direction Method of Multipliers (ADMM) algorithm. The ADMM 
method is known to be particularly slow, see [3|. A method of acceleration based 
on the Nesterov method, [9], is described in [5| and we propose a modification of 
this to increase the speed of the algorithm from m- The algorithm is a variant of 
the gradient descent method with an over-relaxation step. 

Let / and g be functions and that we wish to minimize f{u) -\- g{v) subject to 
Au -1- Bv = c. Let p be a penalty parameter as in [3] . The accelerated algorithm is 
given in Algorithmic To compute the compressed manifold modes we proceed as 
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Algorithm 1 Accelerated Calculation of CMMs. 
Require: oi = 1, = 0.999. 

1: for /c = 1, 2, 3,... do 


2: 

Uk 

= arg min 

f(.u) + ^\\Au + 

Bvk + c- 

-vii^ 

3: 

Vk 

= arg min 

9{v) + ^\\Auk 4- 

- Bv + c - 

■ "^fe P 

4: 

Afe 

— Xk E Auk + Bv]^ + c 



5: 

Cfe 

II 

-Xk\\^ + \\B{vk 


6: 

if « 

Cfe < 

then 



7: Q^/c = 1 

8: end if 


9: 

10 : 

11 : 


1 + \J\ + 4q:| 

«fc+i --- 


l^/c+l — l^/c + 
^/c+1 = + 


H/c ~ 1 

H/c+l 
life -1 

Hfe+I 


- Vfe-l) 

(Ajk — Afc-i) 


12: end for 


described in Section 3.2 of m- The optimization function is split into the sum of 
three (rather than two) functions: 

IV($^L$)+M||$||i + t($) 


where the indicator function l is defined by 


r 0, if = Id, 

[ oc, otherwise. 


They then reformulate the problem as 


min + Tt{E'^LE) + /ills'll! such that ^ = S, ^ = E. 


We proceed similarly to m and use the following in Algorithm 1 : f = u = 

■ E ■ 


V = 


We set 


and 


( 

■ E ' 

w 

■ Ti{E'^LE) ' 

1 

_ 5 _ 




A = 


B = 


-I 0 
0 -I 


and c = 0. 


Note that there is no guarantee that this algorithm converges. That the algo¬ 
rithm in [5] converges for a weakly convex optimization problem is shown in [5]. 
The essential part of the proof is the monotonic decrease in a certain residual. The 
optimization problem we wish to solve is not weakly convex and the residual does 
not monotonically decrease (as can be seen in examples). Hence we make a modifi¬ 
cation to the restart rule to ensure that the algorithm does not become stuck on the 
same values. The acceleration process requires a parameter, denoted r^, to initiate 
the restart. In all the experiments in this paper r] was set to 0.999. 

In the implementation used for the experiments in this paper (and in [TO]) the 
penalty parameter is varied as in Section 3.4 of [3]. This greatly increases the speed 
of both algorithms. 
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Furthermore, convergence was determined by the use of the primal and dual 
residuals. The primal residual at iteration k is = Au^ + Bv^ — c and the dual 
residual is Sk = pA^B{vk — Vk-i). The precise condition for stopping the algorithm 
are described in Section 3.3.1 of [3]. The algorithm terminates when ||r/c ||2 < 
and ||5 /c|| 2 < for the two tolerances, and These are determined from 

eP" = ^ ^rei max 1 11$| I 2 , ^\\E\\l + \\S\\l^ , 

gdual ^ ^gabs_^grel||;^||^ 

where n is the number of vertices for the mesh, is an absolute tolerance, e*'®' is 
a relative tolerance and A is the dual variable in the algorithm. In all experiments, 
following [To], we took = 10“^ and = 10“^. 

In the experiments the two algorithms were given the same random initialization 
of numbers between 0 and 1. This was repeated 30 times for each mesh. Table 
compares the average number of iterations and average time taken for the acceler¬ 
ated algorithm described above and the one from m- The meshes Stand, Crouch 
and Bent Limbs are the those from first, second and third lines of Figure [^respec¬ 
tively; L-shape is the L-shape from Figure Aquarius is the water carrier statue 
mesh from Figure The experiments were performed on an iMac with 3.4 GHz 
Intel core i7 and 8GB RAM. 

In the experiments the accelerated algorithm required on average 47% fewer 
iterations than the original algorithm with a range between 9% and 88%. The im¬ 
provement in speed was nearly as good with an average 44% improvement. The 
range of improvements was 2% to 87%. Hence, we can conclude that the improve¬ 
ments are significant and the algorithm is worth implementing when compared to 
the original. 

Typical graphs comparing the sums of the prime and dual residuals for the two 
algorithms are shown in Figure These were selected to give some insight into how 
the accelerated version behaves. We can see that its residual sum is a condensed 
and bumpier version of the unaccelerated version’s. The combined residual for the 
accelerated algorithm looks like it has greater fluctuations and this is indeed the 
case. This is pictured in close up Figure where one can see that these are Nesterov 
‘ripples’ - a common feature of Nesterov methods. The algorithm restarts at the 
top of the bounce. 

4 A natural order for compressed manifold modes 

For the Laplace-Beltrami operator we can naturally order the eigenfunctions by their 
associated eigenvalues. We shall now describe the correct analogous natural ordering 
for compressed manifold modes. One could naively order the modes according to 
their Dirichlet energy. That is, ignore the additional Li term in the minimisation. 
This leads to poor results as one can see by comparing Figure (our ordering) 
and Figure]^ (Dirichlet ordering). The former is better than the latter for these 
near-isometric surfaces. 

In Figure the modes for a mesh fall naturally into pairs of consecutive modes 
with the exception of the third row where one could argue that positions 4 and 

5 should be exchanged. One could also argue that the head- and torso-supported 
modes are not in the same sequence in the three meshes. However, they do only 
occur in positions 7 and 8. Contrast these small differences with the larger ones in 
Figurethe Dirichlet ordering. Here the head-supported mode appears in positions 
3, 4, and 8. The torso supported mode appears in positions 5 and 6. Furthermore, 
no mesh consistently has what we could call natural pairs and in particular the 
bottom mesh has most of the modes not occurring in consecutive pairs. 
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Stand, /i = 0.0075, K = Q Stand, /i = 0.001, K = Q 




Aquarius, fi = 0.0075, K = 6 


Stand, fi = 0.0075, K = 10 



Crouch, /i = 0.0075, K = 10 


Figure 3: Comparison of typical sum of prime and dual residual for accelerated 
ADMM (blue) and ADMM (red). 
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Mesh 


K 

Fast 

ADMM 

iter’s 

ADMM 

iter’s 

Percent 

reduc¬ 

tion 

Fast 

ADMM 

time 

(in 

sec) 

ADMM 

time 

(in 

sec) 

Percent 

reduc¬ 

tion 

Stand 

0.008 

10 

1381 

2889 

48 

43 

86 

45 

Crouch 

0.008 

10 

1140 

3177 

57 

38 

89 

49 

Bent limbs 

0.008 

10 

1218 

2373 

49 

38 

71 

47 

Stand 

0.001 

6 

1329 

1980 

33 

37 

53 

30 

Crouch 

0.001 

6 

1916 

2704 

30 

48 

66 

27 

Bent limbs 

0.001 

6 

1813 

3640 

46 

50 

95 

43 

Stand 

0.0075 

6 

829 

1498 

44 

23 

40 

42 

Crouch 

0.0075 

6 

902 

1660 

44 

23 

40 

41 

Bent limbs 

0.0075 

6 

922 

1645 

43 

25 

44 

41 

L Shape 

0.02 

10 

7867 

18597 

57 

229 

469 

51 

Aquarius 

0.0075 

6 

1568 

4349 

60 

243 

661 

59 

Aquarius 

0.0075 

10 

2220 

5182 

55 

412 

875 

50 

Average 





47 



44 


Table 1: Comparison of the accelerated and original ADMM algorithms. For each 
mesh the average of 30 iterations is taken. 



Figure 4: A close up (with scale) of the Nesterov ripples. Mesh was Stand with 
11 = 0.008 and K = 6. 






























Having shown the superiority of the new ordering, let us now describe it in 
detail. Let L = A~^W denote the Laplace-Beltrami operator where A is the mass 
matrix and W the weight matrix. Suppose we wish to find K modes. Then the 
modes can be represented as the K columns of ^ where ^ is determined by the 
constrained optimization problem 

nnnTr + /i| |^| |i such that A^ = Id . 

Since A is symmetric the constraint A^ = Id determines K{K — 1)/2 distinct 
equations. If we apply the Lagrange multiplier method we require K{K — l)/2 
Lagrange multipliers with 1 < i < k and i < j < k. 

Lemma 4.1 The Lagmngian function, C, for this constrained optimization problem 
is 

£($) =Tr(i>'^l^i>) +/u||$||i -Tr((i>'^yli>-Id) A) 

where A is the real symmetric K x K matrix with entries [A]ij = Xij for 1 < i < k 
and i < j < k. 

Proof. We can define the Hadamard product of matrices X and Y, denoted X oY, 
to be their elementwise product. That is, 


[XoY]ij = [X]i,-[Y]ij, 


where [X]ij is the {i,j) entry of the matrix X. Now, the Lagrangian function is 


£($) =Tr + mII^IIi - Y - Id) o A]... 


hJ 


From Lemma 5.1.5 of we have that ^ [X o Y].. = Tr(XT^). From this the 


result follows. 


□ 


dC 


The method of Lagrangian multipliers requires the solution of — = 0. Using the 
facts that for matrices X, Y, Z and function /, 

f-T t{XYX'^Z) = ZXY + Z'^XY'^, ([T3] equation 118), 

uX 

f) f d 

fix) = ( ^nx) 


dXT- 


dL 


we see that the equation = 0 gives 


2IU4> + /isign(<I>) — 2H4>A = 0 

+ |$'^sign($) = $^A$A 

q>Twq> + 1$^ sign($) = A. 

By considering the diagonal of both sides, and noting that for a vector ip we have 
sign((/:?) = ||(^||i, we can associate a number to each column of 4>, i.e., to a 
compressed manifold mode. 

Definition 4.2 Let ip he a compressed manifold mode of the operator A~^W with 
compression factor fi. The compressed eigenvalue, \, associated to ip is the 
number 

a = £^W"£ + |||£||i. 


9 





Figure 5: Ordering for three meshes, cf. Figure 13 of m- 



Stand 

Crouch 

Bent limbs 

Teapot 

lO-^x 

1 

19.6155 

20.3833 

20.5455 

34.1979 

2 

19.8722 

20.5321 

20.6391 

57.9528 

3 

21.8682 

21.4825 

24.6300 

60.1676 

4 

23.7901 

23.0605 

26.7960 

60.1705 

5 

27.3304 

25.6135 

27.9656 

60.9936 

6 

28.1828 

27.9401 

29.0192 

62.9510 

7 

28.7432 

30.4192 

29.0294 

63.0666 

8 

34.6949 

33.1840 

33.6681 

63.2446 

9 

39.4302 

36.5946 

38.3945 

63.3030 

10 

39.5416 

36.8892 

38.4758 

63.3360 


Table 2: Compressed eigenvalues for some meshes. 


Note that when ja = 0 and K is the the number of rows of W we get the standard 
eigenvalues of the operator A~^W. 

Figure 13 of [10] demonstrated that their algorithm consistently found similar 
modes for non-isometrically deformed meshes, in particular for a standing, crouching 
and limb bending man from the SCAPE dataset. This experiment was rerun with 
the accelerated algorithm (with /i = 0.008) and the modes were ordered and flipped 
(see Section]^ for an explanation of the latter). This is pictured in Figure]^ and the 
resulting compressed eigenvalues are given in Table 

The compressed eigenvalues can be closely grouped. For example, Aquarius, 
the Water Carrier, pictured in Figurehas compressed eigenvalues 2.2223, 2.4737, 
2.5792, 2.6721, 2.8335, and 3.4063 for = 6 and /i = 0.001. (The accelerated 
algorithm took 3308 iterations from a random initialization.) 

For the human meshes we can see in Table that they range from 19.6155 to 
39.5416. To put this into perspective, the first ten eigenvalues of the LBO range 
from 0 to 31.7886. One can see clearly in Figure that eight of the modes occur 
in obvious pairs: pair of feet, pair of legs, pair of hands, and pair of arms. This 
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Figure 6: Ordering for three meshes by only the Dirichlet energy, note that the 
ordering is not as good as Figure 


pairing is reflected in the compressed eigenvalues, for example for the Stand mesh 
the compressed modes supported on the feet correspond to compressed eigenvalues 
19.6155 and 19.8722. This means that it is highly likely that, just as in the eigen¬ 
value case of symmetric objects, numerical errors can cause compressed eigenvalues 
to be out of order. This is possibly what is happening with the human poses where 
the first compressed eigenvalue is sometimes the left foot, sometimes the right. 

A set of twenty modes is shown in Figure for the classic teapot. Here fi = 
0.0008 and the initialisation was not random but the first twenty eigenfunctions of 
the discrete Laplace-Beltrami operator. Note again that the supports of the modes 
coincide with what a human might identify: the teapot spout, handle, lids, and knob 
of lid. The modes divide the spout in three. It would interesting to investigate how 
varying the parameter fi affects symmetry in the support of the modes. 

Note again, that similar modes arising from symmetries have closely grouped 
compressed eigenvalues. This can be seen in Tablewhere the 3rd to 5th compressed 
eigenvalues correspond to the three regions on the teapot lid. 


5 Accuracy and Consistency 

A measure of the accuracy of a numerical approximation of an eigenvector u of a 
linear map L can be calculated by Lv — Xv. The theory behind the calculation of 
ordering gives us an analogous measure of the accuracy of the algorithm in calcu¬ 
lating modes. A mode cp with compressed eigenvalue should satisfy the condition 
W(p ^ {fi/2) sign((/:?) = XAcp and this can be used to check accuracy. 

In both the original and accelerated versions of the algorithm the average error 
of entries of the vector Wcp (/i/2) sign((/?) — XAcp was of the order 10“^ to 10“^. 
Hence, although we cannot have much confidence in the accuracy of the methods 
we can say that accuracy is not lost by using the faster algorithm. The accuracy 
aspect of CMMs requires further investigation. In particular, what is the trade-off 
between speed and accuracy and what is the relationship between the algorithm 
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Figure 7: Modes on the teapot. 


stopping tolerances and and accuracy. 

One method to improve the accuracy was to perform 200-400 iterations of the 
accelerated algorithm, set all entries in 4> below a preset tolerance to zero and 
then restart the algorithm with this new 4>. This did not conclusively improve 
accuracy but had the advantage of occasionally increasing the speed significantly. 
The results were inconsistent and further study is required. One could also try a 
similar approach by setting to zero the entries of 4> that correspond to large absolute 
entries of Wcp {ja/2) sign((^) — XAcp for the different (p. 

We now turn to consistency. Given a random initialization it is clearly plausible 
that resulting modes are inconsistent, that is distinct runs of the algorithm with 
fixed /i and K, etc, may produce a different collection of modes. However, experi¬ 
ments show that for certain values of p the calculation of modes by the accelerated 
method is very consistent. For the SCAPE dataset models (Stand, Crouch and 
Bent Limbs in this paper) repeated experiments on random initializations involving 
numbers from 0 to 1 produced less than 1% inconsistent results for fi around 0.008. 
For other models, for example Aquarius, it was harder to locate a good value of p 
to use. 

The value of p is important for applications as p determines the ‘size’ of the 
support of modes, that is the amount the support covers the model. If p is large, 
then the supported area is small. If the supported area is small, then it is likely that 
supported areas do not interact during the iterations and this can in theory lead 
to inconsistency. Hence, the value of p is important for consistency but is poorly 
understood. 

One way of avoiding inconsistency is to use the eigenfunctions of the Laplace- 
Beltrami operator as the initialization of 4>. This need not give the same K modes 
as the random initialization and it would be interesting to investigate further the 
use of this initialization as it worked quite well in some experiments. For example, 
the teapot in Figure was calculated using the first twenty eigenfunctions and this 
was very consistent for small changes of p in the sense we had the same modes in 
the same order. Nonetheless, what are the quantitative changes to the modes when 
fi is varied is still an open question. 

When one changes the number K of modes to be calculated, then one gets 
inconsistent results. For example, compare Figure with Figure where K is 
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Figure 8: Effect of the choice of K. 


8 and 10 respectively. The modes calculated for = 8 are not the first 8 modes 
calculated for K = 10. However, Figurej^gives one confidence that good consistency 
is achievable even between near isometric shapes. 

6 Orientation and flipping of modes 

A disadvantage of the Laplacian unit eigenfunctions is that they are defined only 
up to sign and there is no obvious way to ‘flip’ them so that eigenfunctions between 
different manifolds can be directly compared. For these eigenfunctions one expects 
an ‘equal amount’ of positive and negative values. That is, the integral of the 
positive values equals the integral of the negative values. (This is because the 
eigenfunctions are orthogonal to the constant vector, the eigenfunction for the zero 
eigenvalue.) However, for general compressed modes there is not this balancing of 
the positive and negative. Hence we can flip modes so that the ambiguity of parity 
is removed.. 

In practice, the values of a mode are dominated by numbers of either positive or 
negative values. For a mode ip the number sign {mdix{ip) + min((/9)) will, in practice, 
be ±1 depending on the sign of the dominant values. One can then use this to 
change the sign of the columns of In this way, the modes chosen to be dominated 
by positive values. A problem only occurs when max((/:?) = — min((/9) but this is 
unlikely in practice. 

An alternative and perhaps more theoretically sound way to flip the mode is 
to integrate it over the manifold and define its sign to be the sign of the resulting 
value. 

That the positive values dominate after flipping can be seen in the figures where 
red is positive. 

7 Lumped and unlumped mass matrices 

The mass matrix Am L = A~^W contains information about the area around the 
vertices. Many variations are possible, for example one-third area, Voronoi area or 
uniform area, see [8]. Given such a matrix one can ‘lump’ the matrix by summing 
all the elements of a row and putting the result in the diagonal. A mass matrix 
with non-zero entries only on the diagonal is called a lumped matrix. 

The algorithm in HO] is given for lumped matrices. In that paper D is a diagonal 
matrix and it is easy to calculate and its inverse where is formed 

from the square roots of the diagonal entries. For a positive definite matrix B such 
as the mass matrix there exists a square root matrix such that B^I^B^I^ = B. 
However, the calculation of this can be computationally expensive and hence though 
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it is theoretically possible with this method to include unlumped mass matrices in 
the algorithm of m this may seriously impact the speed of computation. 

Nonetheless, unlumped matrices can be used. In [TO] the first part of the algo¬ 
rithm involves letting Y = — Us E Ue)^ finding an SVD factorization 

of Y^Y = VWV^ and then the closed form of the minimization problem is 


If instead we tak^F =_S-UsEE^ Ue, and so F = D^/‘^Y then F^F = VWV^ 
simply becomes, Y^DY = VWV^ and the closed form of the minimization problem 
is then _ _ 

$ = £>-1/2 ^£)l/2y^ y^^\|2yT ^ YVW'^I'^V'^. 

Thus with this modification to the algorithm we can avoid the expensive computa¬ 
tion of the square root of the mass matrix and allow use of unlumped matrices. 

Numerous experiments with different meshes, values of K and /i, failed to show 
that either the lumped or unlumped version was superior in terms of speed or 
accuracy. It seems unlikely that lumping does not at least have some effect and 
so while the experiments were inconclusive it may be the case that there do exist 
identifiable situations in which one method is superior to the other. 

To ensure consistency and comparability the experiments described in this paper 
were done with lumped matrices. It should be noted that lumped matrices are used 
for the LBO much more commonly than unlumped. The aim here is to show that 
we can use unlumped if we wish. 

8 Conclusions and future work 

Compressed manifold modes were demonstrated in HO] to be robust with respect to 
noise and holes and have the potential for use in geometry processing applications. 
Indeed as they generalize to surfaces the compressed modes introduced in El they 
also have potential for application in solving PDEs on surfaces. In this paper it has 
been shown that the algorithm of m for the computation of compressed manifold 
modes can be considerably improved - by almost a factor of 2. The correct method 
for ordering modes has been demonstrated (with proof) and the modes have been 
oriented. 

The study of the numerical calculation of eigenvalues has had the advantage 
of decades of work and the development of many optimization techniques and so 
unsurprisingly the calculation of compressed manifold modes is not competitive in 
terms of speed. This paper helps reduce this uncompetitiveness but the development 
of further improvements would be welcome. 

Compressed manifold modes are of interest because they are functions with local 
support. A crucial insight is that if one is looking for a collection of functions to 
form a basis of functions on the manifold, then speed and accuracy are not that 
important, it is the ortho normality of the collection that is important. In the case 
of the CMM algorithms the locality of support and orthonormality can, and in 
the experiments usually do, appear rapidly, say within 400-600 iterations of the 
algorithm. Therefore it would be good to analyse how quickly these properties are 
achieved. 

Further work of interest includes the following: 

(i). Investigate and define for numerical calculations what it means for the modes 
to be locally supported and orthonormal. For the latter we need to check that 
is close to the identity. How close for practical applications? 
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(ii) . The random initialization can lead to a wide variation in computation time. 

For example, in the 30 calculations used in Table the Bent Limbs mesh 
with K = 10, jii = 0.008 took between 736 and 4323 iterations, the latter is 
a factor of nearly 6 greater than the former. Can we do better or at least be 
more consistent? 

(iii) . Can inconsistency of the collection of modes calculated be reduced by restrict¬ 

ing the band of random numbers used as an initialization? See Section 

(iv) . Can accuracy of the algorithms be improved by some method, for example 

truncating the support as in Section 

(v) . How do changes in the compression factor /i affect the consistency and accu¬ 

racy of the algorithms? 
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